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EXECUTIVE SUMMARY 


The general aim of this Phase I SBIR project was to extend and refine 
existing computer simulation capabilities in order to compute equilibrium and 
nonequilibrium properties of concentrated colloidal suspensions. The immed¬ 
iate application of this work is to guide improvements in the characteristics 
of nonaqueous colloidal fuel slurries that are currently being developed by 
the Air Force for use as propellants. 

Equilibrium and nonequilibrium Brownian dynamics simulation techniques 
were employed in this work. The simulations include the effects of inter¬ 
particle forces, Brownian motion and shear flow on the motion of particles. 
Hydrodynamic interactions were neglected in accordance with the Phase I work 
plan. 


Equilibrium and nonequilibrium Brownian dynamics (NEBD) simulations were 
performed for concentrated aqueous and nonaqueous colloidal suspensions. 

These are the first NEBD simulations in which the shear viscosities of suspen¬ 
sions have been calculated. Stability estimates were also made for sterically 
stabilized nonaqueous suspensions. These results can be understood in terms 
of the changes occurring in the repulsive steric potential and attractive 
van der Waals potential as the adsorbed polymer concentration and Hamaker con¬ 
stant are changed. 

These calculations provide quantitative evidence for the disturbance of 
the equilibrium structure of the dispersions by shearing. Shear viscosities 
of suspensions have been calculated as a function of shear rate and particle 
volume fraction. Good agreement was obtained with experimental viscosities 
for comparable systems. The simulation results are sensitive to the type, 
strength and range of the potential interactions used to describe the disper¬ 
sion. Thee overall conclusion of this Phase I study is that the NEBD tech¬ 
nique is capable of describing the essential features of sheared suspension 
behavior. These Phase I results establish the feasibility and desirability of 
using NEBD simulation methods as a predictive engineering tool for the design 
of slurries. 
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1. INTRODUCTION 


The general aim of this Phase I SBIR project was to extend and refine 
existing computer simulation capabilities in order to compute equilibrium and 
nonequilibrium properties of concentrated colloidal suspensions. Nonaqueous 
colloidal fuel slurries are currently being developed by the Air Force for use 
as propellants. The primary goal of this research was to perform, for the 
first time, nonequilibrium Brownian dynamics simulations of sheared concen¬ 
trated suspensions in order to calculate the shear viscosity of these systems. 
This has been accomplished. By comparing the calculated viscosities and their 
dependence on particle volume fraction and shear rate with observed properties 
of real systems, the feasibility of this approach as a slurry design tool has 
been established. Sample calculations of suspension stability using classical 
coagulation theory have also been done to illustrate the interplay of repul¬ 
sive and attractive forces in stabilizing suspensions. 

Because of the many different types of forces and interactions between 
colloidal particles, these suspensions display complicated behavior, and 
current theoretical understanding of them is limited. 1-5 Exact theories do 
not exist except for a few idealized, limiting cases. Approximate theories 
are invariably forced to omit or oversimplify one or more of the important 
classes of particle interactions for the sake of mathematical tractability. 
Computer simulation techniques, on the other hand, are capable of treating 
more-or-less exactly all forms of particle interactions, at the expense, 
admittedly, of ever increasing computational effort with growth in the size of 
the system or number of effects treated. However, modern computers are 
sufficiently fast to overcome these limitations for many systems of practical 
interest. Thus, these methods are a promising means for obtaining insight 
into the microscopic reasons underlying the observed behavior of concentrated 
suspensions. 

Equilibrium and nonequilibrium simulations of concentrated aqueous and 
nonaqueous suspensions were performed. Aqueous systems were studied for two 
reasons. First, results of prior equilibrium simulations of aqueous systems 
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have been published and are available for comparison. Second, experimental 
shear viscosity data on well-characterized suspensions cover a somewhat 
broader range of conditions for aqueous systems, which have been studied more 
extensively than nonaqueous systems. Equilibrium simulations of strongly 
interacting dilute systems were also performed under previous PSI in-house 
support. The results of these comparisons are described in Appendix A as 
additional background material. 

It should be noted that these are the first nonequilibrium Brownian 
dynamics (NEBD) simulations in which the shear viscosity of a suspension has 
been calculated. Several different interparticle force laws were studied, 
although hydrodynamic interactions were neglected in accordance with the lim¬ 
ited Phase I work plan. Previous BD simulations with shear and hydrodynamic 
interactions,® have been concerned with the dynamics of coagulation and 
deflocculation of isolated pairs of particles. Nonequilibrium Stokesian 
dynamics simulations were developed by Bossis and Brady 7 to treat suspensions 
of particles large enough that Brownian motion of the particles can legiti¬ 
mately be neglected. This method includes hydrodynamic interactions but has 
so far only been used to calculate the shear viscosity of two dimensional sus¬ 
pensions. 7 ’ 8 Many nonequilibrium molecular dynamics (NEMD) simulations have 
been performed to compute the shear viscosity of molecular fluids. 9-12 
Because of the vast differences in the size and time scales governing molec¬ 
ular and colloidal particle motion, the NEMD simulations must be performed at 
extremely high shear rates, exceeding by many orders of magnitude values that 
can be attained experimentally. Although NEMD techniques are an extremely 
important development in computational statistical mechanics, they are incap¬ 
able of treating many inherently colloidal phenomena such as viscous damping 
of particle motion, hydrodynamic interaction between particles or Brownian 
motion. Thus, NEMD simulations can provide only some qualitative insights 
into the rheology of particulate suspensions. 

Before describing the results of the nonequilibrium simulations, the 
theoretical background underlying the computational method is presented in 
Section 2. Following this, the interparticle pair potentials used in the 
simulations are discussed in Section 3. As noted, simulations were performed 
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with pair potentials appropriate for aqueous as well as nonaqueous disper¬ 
sions. The results of the equilibrium and nonequilibrium simulations of con¬ 
centrated dispersions are then presented and discussed in Section 4. The 
stability ratio for nonaqueous sytems is also discussed. Appendix A contains 
the results of the equilibrium simulations of dilute, but strongly interacting 
aqueous dispersions. 
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2. THEORETICAL BACKGROUND 


2.1 GENERAL CONSIDERATIONS 

The Brownian dynamics (BD) simulation technique is employed in this work. 
This method is an outgrowth of the molecular dynamics (MD) simulation tech¬ 
niques that are now a well established part of the methodology of modern 
statistical mechanics. 13 The MD method involves numerically solving Newton's 
equations of motion for all N particles in the system. It is a deterministic, 
rather than a stochastic, technique. By keeping track of all the particles' 
trajectories, a representative sampling of phase space is generated, and aver¬ 
ages of system properties can be calculated. Both dynamic and static system 
properties can be treated. The method is limited to relatively small numbers 
of particles, generally less than a thousand, because of computer speed and 
size limitations. 

Because of the vast difference in size and mass between a typical suspen¬ 
ded particle and a molecule of the suspending fluid, the time and length 
scales for simultaneous motion of molecules and particles cover too many 
orders of magnitude to be directly treated by MD techniques. To overcome this 
problem, BD techniques 14-18 have been developed to simulate suspensions. In 
the BD technique, the suspending fluid is treated as a hydrodynamic and die¬ 
lectric continuum and the motions of the suspended particles are described by 
Langevin equations. Langevin equations are stochastic differential equations 
that play the role of Newton's equations in the molecular dynamics simulation. 
The Langevin equations describe the motions of the particles subject to the 
action of interparticle and external forces, frictional forces arising from 
solvent drag, and random Brownian motion forces. The Brownian forces are ass¬ 
umed to act on a time scale that is much shorter than that required to produce 
even a modest change in the particle's location or velocity, and they must be 
treated statistically, as described below, to account for their cumulative 
efffect on the motion of the colloidal particle. From the solution of the 
Langevin equations for the interacting many particle system, positions of the 
N colloidal particles can be computed as a function of time. From this 




information, one can calculate fundamental quantities of interest such as mean 
square particle displacements, pair distributuion functions, and transport 
properties such as the self-diffusion coefficient and the shear viscosity. 


2.2 PARTICLE INTERACTIONS 


Solution of the equations of motion for individual particle displacements 
requires knowledge of the interparticle forces and interactions affecting each 
particle. These may be broken down into three groups. 

First, there are random Langevin forces responsible for the Brownian 
motion of the particles. Langevin forces arise from the (small) fluctuations 
in the (very large) number of collisions experienced by different areas of the 
large particle with the much smaller molecules of the suspending fluid. 

The second class of forces consists of the "direct" interparticle forces. 
These include: (1) repulsive (or excluded volume) forces due to the finite 
size and relatively small deformability of the particles, (2) short-ranged 
oscillatory forces between closely spaced particles due to fluid structure and 
packing, 1 ^* 20 (3) forces due to absorbed polymer molecules which, under con¬ 
ditions giving rise to net repulsion, result in so-called steric stabili¬ 
zation, 21 (4) attractive London-van der Waals forces, and (5) repulsive forces 
due to the overlap of electrical double layers surrounding similar particles. 
The latter two forces form the basis for the classical DLVO theory of colloid 
stability. 22 * 22 The attractive van der Waals force between large suspended 
particles at small separations d varies as d -2 , 24 considerably different from 
the short ranged d -7 dependence typical of molecular interactions. Electrical 
double layer repulsions (or attractions for dissimilar particles under the 
right conditions) result from overlap of the diffuse clouds of counter ions 
surrounding a charged suspension particle. For situations of weak overlap the 
interaction force for large particles varies exponentially with d; with the 
length scale set by the familar Debye parameter 25 (or inverse screening 
length), <. In non-aqueous media, electrical double layer interactions are 
generally unimportant. 
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The third group of interactiona consist* of hydrodynamic interactions 
between particles. The term refers to the interaction between particles due 
to the velocity fields set up in the fluid by the moving particles. These 
velocity disturbances propagate through the fluid, perturbing the motion of 
nearby particles. 

Hydrodynamic interactions ars neglected in this early developmental work, 
and will only be cursorily treated in describing the algorithm for colloids, 
particle motion. 

2.3 LANGEV1N EQUATIONS 


The Langevin equation governing the motion of colloidal particle 1 ot 
mass m^ can be written 6 as 


d ?i 
“i dt~ 


- i 
i ~ ij 


(v - U ♦ 




♦ r ♦ l 

-i "i 


where v< is the particle velocity, Ft is the total force ldire< . 
plus external, if present) on the particle, is tha random U ngevm totem 
responsible for the particle's Brownian motion, u^ is tha undisturbed tin* 
velocity at the location of particle 1 due to an externally imposed tlow 
field. Hydrodynamic interactions batween particles are governed ay the con¬ 
figuration dependent second rank friction (or resistance) tensor ^ and tha 
configuration dependent third rank shear tensor S. The rate of strain tensor 
E is given as usual by the expression 


£ “ i&i ♦ (Zh) T 1 • < 

For the case of simple shear flow in the x direction, the unperturbed fluid 
velocity flow field is given by 

u, ■ v y, e ( 

~i ' i ~x 

where y i* the shear rate (s -1 ) and y^ is tha y coordinate of particle i in 

some Cartesian coordinate system with unit basis vectors (e x , e y , e s ). 
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The first term in Eq. (1) rsprsssnts ths frictional resistanca on tha 
particla due to aolvant drag, hydrodynamic intaractions and convactiva affacts 
in a ahaarad auspansion. With tha naglact of hydrodynamic intaractions, tha 


shaar tanaor vaniahaa, 


2i 


(4) 


and only tha diagonal coaponanta of tha friction tanaor raaain, 




‘i 5 i d i 


(S> 


Hara, 1^, tha single particla friction coefficient, is customarily evaluated 
using Stokes' law, 


; i 


San a 


( 6 ) 


where is the pure fluid viscosity (Pa*s) and a is the radius of the parti¬ 
cle. Because of Eqs. (4) and (5), with the neglect of hydrodynamic interac¬ 
tions tha Langevin Equations for particla motion simplify considerably to the 
form 


*i dt 


M~i 


F ♦ L 
~i ~i 


(7) 


The complicated configuration dependence of F^ still precludes an exact solu¬ 
tion of these equations. In order to proceed with the numerical simulation of 
particla motion, an algorithm is required for computing the change in position 
of each particle during a small time step at. Heuristic arguments will be 
used to obtain algorithms without and with shear flow that are special cases 
of those obtained with more rigor by Ermak and HcCammon 17 and Ansell, 
Dickinson, and Ludvigsen. 6 


Particles undergoing Brownian motion incur significant changes in posi¬ 
tion only on a time scale that is long compared to tha typical time for tha 
particle's instantaneous velocity to become uncorrslated with its value at any 
earlier time. 26 Velocity fluctuations are rapidly damped by fluid friction, 
and all of the forces acting on a particle are nearly always in balance. The 
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velocity relaxation time ecele is set by the ratio m^/^. Thus for a tine 
interval At >> m^/c, i. the left hand aide of Eq. (7) may be neglected, and the 
remaining terms may be solved for v A to give 

»i * * £ l /; i ♦ i t /c t • (8) 

Equation (8) is the i'th member of a set of N first order, coupled 
stochastic differential equations for the poaitions r^ of the N colloidal 
Brownian particles in the suspension. Because dr^/dt • v 1# it is apparently a 
simple matter to solve Eq. (8) for the displacement Ar incurred in a time 
interval short enough that the total force F^ acting on the particle remains 
constant. However, solving a stochastic differential equation is different 
from solving an ordinary partial differential equation. Given the stipulation 
that At >> m^/;^, the product L^At has no meaning because varies extremely 
rapidly and randomly within the interval At. A correct integration of Eq. (8) 
results in the following prescription for Ar^: 

Ar - [u ♦ D F Al)At ♦ R. (9) 

-l '-'1 1 —1 

where the random displacement experienced by the particle in At is defined as 

t*At 

5 i - ; t/^] / L i dr (10) 

and the self-diffusion coefficient, 0^, of a particle in a highly dilute 
suspension is given by the Stokes-Einstein expression 

Di - kT/^i - kT/( 6 wn 0 ■> (11) 

In order to complete the algorithm, the statistical properties of Rj must be 
specified. These follow from two basic assumptions of Brownian motion theory 
r * 9 * r<J ing the behavior of L*. First, the average of over all possible 
fluctuations is sero, 

<L > - 0 . Hi) 

■*1 L ■'* 
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Second, values of at different times are uncorrelated, i.e., behaves as 
a white noise source, 


‘ii**!* 


W'l 




5 (t. 


V 


( 13) 


The subscript L on the angle brackets in Eq. (12) and (13) indicates an 
average over all the fluctuating values of L^. The constants preceding the 
delta function in Eq. (13) are a consequence of a form of the Fluctuation- 
Dissipation Theorem. They are there to insure that the long time average 
value of v^ 2 , as determined by exactly solving Eq. (7) without shear (u^-0), 
equals 3kT/m^, the equipartition value of equilibrium statistical mechanics. 


Using the properties shown in Eq. (12) and (13), it is easy to see that 


<R. > - 0 

-i L 


( 14) 


and that 


<R > - 6D.it . 

1 L X 


(IS) 


Equation (9) may now be written as 


1/2 


A ~i " ^“i * D i£i /kT ) At ^ 


( 16 ) 


where the components of Q< are Gaussian random numbers whose mean and mean 
square values are zero and one, respectively: 


Qxi * Q xi " 1 * (17) 

Equations (16) and (17) comprise the algorithm for computing the dis¬ 
placements of the colloidal particles. The displacement is seen to result 
from the superposition of a term linear in At, which is due to the shear field 
(if any) and direct particle interactions, with the Brownian jump term, which 
varies as the square root of At. These equations also give the correct expres¬ 
sions to first order in At for the mean and mean square displacements of a 
diffusing particle in a shear field and subject to interparticle or external 
forces (but without hydrodynamic interactions): 


9 




<A* > - (u A ♦ D^F^/kT) At 


( 18) 


< (Ar i ) 2> L - 6D^t . (19) 

The general procedure for carrying out a simulation is as follows. Given 
a set of particle coordinates at a time t, the forces on each particle are 
determined, and a set of random numbers Q< is generated. The particle posi- 
tiona are then advanced, and the entire process is repeated. The questions of 
boundary conditions and how to correctly determine the total force on a 
particle will be dealt with in the following sections. 

2.4 PERIODIC BOUNDARY CONDITIONS 


In order to insure that a system of N particles (where N is usually less 
than several hundred) in a cubic cell of volume V (adjusted to give the 
desired density or volume fraction) will represent a bulk system, the primary 
cell is surrounded by 26 image cells each containing N particles in the same 
configuration as the primary cell. Such periodic boundary conditions allow 
particles to pass through cell walls while preserving constant particle number 
(and density) in each cell. 

The aide length l of the cube is determined by specifying the volume 
fraction $ for the suspension in terms of the particle number density and 
volume, 


* - (N/f 3 )(*d 3 /6) (20) 

where d is the particle diameter. 

Consider first the oase without shear flow. Beca tse all of the images of 
a given particle execute the same movements as the particle itself, when the 
particle passes through one of the cell walls, an appropriate image particle 
automatically passes through the opposite wall to replace it in the cell. The 
situation is illustrated in Fig. 1. In practice it is unnecessary to save 
the coordinates of the image particles because these are easily generated by a 
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Figure 1. Two-dimensional illustration of how periodic boundary conditions 
simulate an infinite equilibrium system . Circles denote 
positions before a jump of (exaggerated) magnitude + . Squares 
indicate positions after the jump. Open symbols represent 
images surrounding the primary particles (filled symbol). 

series of reflections, x^+vi, y^+vi, z^^vt, where v - -1 , 0, and +1. For 
example, with the coordinate origin at the center of the box, if particle i 
moves outside the box at x«f/2, i.e., if x i >i/2, it is replaced by an image 
particle now located at (within the cube). Similarly, if then 

x^ would be replaced by x^+i. Similar adjustments are made to y and z 
coordinates whenever necessary. 

In computing the total force on a particle due to interactions with all 
N-1 remaining particles, the minimum image convention^ 7 is used. The cut-off 
distance for particle interactions is taken to be 1/2. Forces between 
particles farther apart than t/2 are set equal to zero. This insures that, 
when rij>t/2, particle i will interact with only the nearest image of particle 
j. The nearest image is located by determining which of the inequalities, 
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A. 







-i/2<x i j<Z/2 

(21a) 

-H/2<y i j<l/2 

(21b) 

-l/2<z i j<l/2 

(21c) 


is violated, followed by adjusting the appropriate components of rj. With a 
few moments of reflection, it should be clear that Eq. (21) will be satisfied 
by every pair of particles, provided the nearest image coordinates of j are 
used whenever appropriate. The situation is illustrated in Fig. 2. 

With a shear flow in the x direction, periodic boundary conditions iden¬ 
tical to those used without shear are imposed for motion in the x and z direc¬ 
tions. Motion in the y direction is subject to the coordinate space version 


O 

O 
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A 

4/ <t/2 

A 
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/ 
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/ r>t/2 


A 

l 

A 

O 

r<l/2^3 

o 

A 

A 

A 


A-5242 

Figure 2. Two-dimensional illustration of how particles interact when 

the minimum image convention is used . The two particles e, A 
in the center square are too far apart to interface. The e 
particle interacts with the nearest triangular image and the 
A particle interacts with the nearest circular image. 
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of the homogeneous shear conditions introduced by Lees and Edwards^ 8 and 
Evans. 10 Image particles in the cubes above the plane y»£/2 and below the 
plane y— 1/2 do not simply replicate the motion in the primary cube. If they 
did, the shear field would exhibit sharp discontinuities at y=±i/2. Because 
of the homogeneous shear condition, image particles lying above y-4/2 or below 
y~-l/2 suffer added x displacements of magnitude 4,yAt in addition to those 
incurred by the primary particles. This insures that image particles with 
larger or smaller y values than those of the primary particles experience 
proportionately larger shear displacements in the appropriate direction. 

The result of this procedure is that after n steps the relative x dis¬ 
placement between a particle and its (original) nearest upper and lower images 
is nfyAt. The particle and these two images lie on a diagonal line in the x-y 
plane that is rotated through an angle a where tana » nyAt. The situation is 
illustrated in Fig. 3. When a particle leaves the box at y«±i,/2, it is 
replaced by a particle with new coordinates (xi±n2.yAt, yi+i); the new x 
coordinate must also be checked and adjusted, if necessary, to insure that it 
satisfies the inequality l/2>Xj i >-l/2. When a equals 45*, the additional x 
displacement between images due to shear vanishes, and the simulation may be 
continued with n reset to zero. 

In computing the total force on particle i, the minimum image convention 
is again used, but allowance must be made for the relative shear displacement 
of image particles. When y^j satisfies the inequality 

-l/2<y i j<l/2 (22) 

the procedure is identical to that in the absence of shear. When y i j>i/2, Xj 
must first be changed to Xj+nfyAt before checking the inequalities in Eq.(21). 
When y i j<-£/2, x^ is first changed to Xj-niyAt before checking Eq.(21). An 
illustration is provided in Fig. 4. 
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Figure 3. Illustration of how homogeneous shear conditions plus periodic 
boundary conditions simulate an infinite system being sheared 
at a shear rate y. Circles mark positions of particles after 
n time steps of At; squares indicate positions after n+1 steps. 
The value of yAt is exaggerated for clarity. The displacement 
of particle • results from the nearly vertical jump due to 
random and interparticle forces plus the horizontal shear 
displacement. After n steps, particle o and its images o lie 
along the solid diagonal rotated through the angle a. After 
n +1 steps, particle ■ and its images □ lie along the dashed 
diagonals that are rotated by the angle 3 . 



Figure 4. Illustration of how particles interact when homogenous shear con¬ 
ditioning and minimum image convention are used. The two particles 
in the center square are too far apart to interact. Each 
interacts with the nearest image of the other, but the effects of 
shear must be accounted for in locating the images. 


jflL_L 


14 




2.5 RADIAL DISTRIBUTION FUNCTION 


The radial distribution function (or pair correlation function) g(r) 
provides information about the local distribution of particles in the fluid 
around any arbitrarily selected particle. Because particles have hard cores 
and interact with each other via other repulsive and attractive forces, the 
local distribution may deviate considerably from the simple bulk average. 
Besides providing insight into the structure of the fluid, g(r) is important 
for another reason. When particle interactions are assumed to be pairwise 
additive, all of the (excess) equilibrium thermodynamic properties of the 
suspension can be obtained solely from g(r). There are a number of simple, 
equivalent ways of defining g(r). If p (*N/V)is the average number density of 
particles in the suspension, pg(r) is the local density at a distance r from 
any particle. Thus, the quantity 4irr 2 pg(r )dr is the number of particles whose 
centers lie within a spherical shell of thickness dr at a distance r from the 
central particle. Finally, N4irr 2 pg(r )dr is the number of pairs of particles 
whose center-to-center distance lies between r and r+dr. Because the 
distances between all pairs of particles are evaluated in determining the 
total force on each particle, this last definition is most convenient in 
numerically calculating g(r) during a simulation. 

2.6 SELF-DIFFUSION COEFFICIENT 

The self-diffusion coefficient for tagged (or tracer) particle motion in 
suspensions is a good indicator of the effects of concentration, interparticle 
forces, and hydrodynamic interactions on the ease with which particles move 
about in the suspension. The time dependent self-diffusion coefficient can be 
computed from the expression 29 * 30 

D - i 4r<Ur(t)) 2 > (23) 

s 6 dt v ' 

where the mean square displacement (m.s.d.) is defined as 
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(24) 


<(Ar(t)) 2 > = |r.(t) - r.(0)| 2 > 

and the angle brackets denote an equilibrium ensemble average. The algorithm 
used to compute the m.s.d. is similar to Van Megen's and Snook's method, 30 
which is based on the operation of a hard-wired correlator used in photon 
correlation spectroscopy experiments. 

The m.s.d. is computed at n s time intervals t of length mAt. The con¬ 
stants m and n s are chosen to give adequate temporal resolution of the m.s.d. 
and to insure that the longest sampling time (n s x) is less than the time at 
which the uncertainty in the calculated value of the m.s.d. becomes too large 
due to statistical fluctuations. Rather than store n s sets of particle coor¬ 
dinates, as Van Megen and Snook do, a suggestion 31 to store the cumulative 
displacements is followed. In this way, one avoids having to correct the 
sequence of coordinates for a particle when that particle passes through any 
wall of the simulation cube. Failure to correct the coordinates would intro¬ 
duce artificial jumps of length l every time a particle moved through a wall. 

The cumulative x, y, and z displacements of the N particles at each of 
the n s times are stored in n s sets of registers. After every m time steps, a 
new set of initial displacements Ar(-r) is generated, ana the total displace- 

»w 

ments corresponding to the largest sampling time n s T are discarded. The new 
displacements are then added to each of the remaining cumulative displace¬ 
ments, thus advancing them by one step of time t. Each of the cumulative 
displacements is then shifted to the next storage register, and the new 
displacements are stored in the first register. In this scheme, the particle 
trajectories are followed backwards through time (which makes no difference in 
calculating the m.s.d.), and a new sequence of cumulative displacements is 
effectively created after every m time steps. From the new sequence of cumu¬ 
lative displacements the squares of the displacements are formed, summed, and 
added to the existing sums for each time. Ultimately, these sums are averaged 
over the total number of starting times for which a sequence of displacements 
was generated. 
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2.7 SHEAR VISCOSITY 


From the nonequilibrium simulations, the shear viscosity coefficient n 
can be obtained by dividing the average shear stress a X y of the suspension by 
the shear rate: 


o m n y (25) 

xy 

With the neglect of hydrodynamic interactions, the average shear stress can be 
calculated by averaging the particle stress tensor:®*^2,33 


a 


xy 



+ 1 
2 1 - 4 > 



i> 3 


y. . f . 

13 X13 


(26) 


The first term is the pure fluid contribution, the second arises from the 
hydrodynamic particle stress in the infinitely dilute limit, and the last term 
represents the contribution of interparticle forces to the stress. Here, the 
angle brackets denote an average over the nonequilibrium particle distribution 
in the sheared suspension. The double sum indicated in Eq. (26) is computed 
after each time step of the simulation. The values are summed and averaged 
over the total number of time steps in the simulation. 


2.8 INITIAL CONDITIONS 


Simulations are generally started with the particles arranged on the 
sites of a fee lattice. Typically, five thousand time steps (or more, if 
needed) are then generated in order to bring the system to equilibrium or to 
steady-state (with shear). Computation of average system properties then com¬ 
mences. In the equilibrium simulations, satisfactory results can usually be 
obtained with five to ten thousand steps. For the nonequilibrium simulations 
much longer runs, of thirty thousand steps or more, are required for averaging 
the particle stress tensor. If the final coordinates from a prior run under 
identical or similar conditions are available, the preliminary startup period 
may be eliminated or significantly reduced. 


3. DIRECT INTERPARTICLE FORCES 


3.1 PAIRWISE ADDITIVITY 


Potential interactions between particles were assumed to be pairwise 
additive. Thus, the total configurational potential energy, U, of N particles 
can be written as 


U 

i 


. 1 . 

i>3 


(27) 


where r^j (= |irj —r^ J) is the distance between the centers of particles i and j. 
The total force on a particle is obtained as usual from the gradient of the 
potential: 


F. = - V.U 

~1 ~i 



z 

j*i 


F. . 
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3.2 AQUEOUS SYSTEMS 

For the aqueous suspensions studied here, only repulsive forces due to 
the overlap of electrical double layers were considered. The range of the 
pair interactions is characterized by the familiar Debye screening length 

..-1 

IS I 


2cN V 2 

■ - •tar) 

where e is the electronic charge (1.602 x 10 -19 C), c is the monovalent 
electrolyte concentration, e is the electrical permittivity of the aqueous 
phase (7.08 x 10” 10 F/m), and N A is Avogadro's number. 
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When the electrolyte concentration and particle diameter are such that <d<3, a 
good approximation for u(rij) is the screened Coulomb (or Yukawa) 
potential,22,29,42a 


u(r) = ire(di|; 0 ) 2 {exp[-<(r-d) ]}/r 


where i}> 0 is the electrical potential at the particle surface. 


When <d>3 a better approximation to u(r) is given by 22 * 29 


u(r)« iredtpo ln{ 1 +exp[ -< (r-d) ] } 


under conditions of constant surface potential. This potential was used in 
simulations of concentrated aqueous dispersions of 1.2 ym diameter spheres 
at two electrolyte concentrations: 10 -4 mol/dm 2 and 10“^ mol/dm 2 . Plots of 
Eq. (30) for these conditions are shown in Fig. 5. 


3.3 NONAQUEOUS SYSTEMS 


For the nonaqueous dispersions studied here, the pair potential was taken 
to be the sum of a repulsive piece arising from steric stabilization and an 
attractive piece due to van der Waals forces: 


u(r) = u s (r) + u A (r) 


Steric stabilization arises from the repulsive interaction of polymer layers 
adsorbed on or bound to particle surfaces. The steeply repulsive potential 
derived by Ottewill and Walker 24 was used: 


u s (r) * 

6 Vp 


(4»~X) C 2R 0 + r )( R o - r )‘ 
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Figure 5. Potentials used for simulating concentrated dispersions . 

where 

R 0 * d + 26 . 

The thickness of the adsorbed polymer layer is denoted by 5, and R Q is 
the distance at which the steric repulsive force is first felt. In the simu¬ 
lation runs, a layer thickness of 0.02 pm was used and a value of 0.013 g/cm^ 
was used for Cp, the concentration of polymer segments per unit volume of the 
adsorbed layer. The remaining parameters were evaluated following Cairns, Van 
Megan and Ottewill-* 5 and are appropriate for poly (12-hydroxystearic acid) 
chains with a dodecane solvent. The density pp of the adsorbed chains is 
1.114 q/cm3 f the difference between the entropic (i|>) and enthalpic (y) inter¬ 
action parameters is 0.3, and the partial molar volume v 1 of the solvent is 
227.3 cm^/mol. 
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A value of 5.0 x 10- 20 J was used for the Haaaker constant, A. A plot of 
Eq. (31) is shown in Fig. 5. 
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RESULTS AND DISCUSSION 


4.1 EQUILIBRIUM SIMULATIONS 

Equilibrium simulations of concsntratad suspensions of 1.2 um diameter 
spheres were performed. Aqueous suspensions at volume fractions of 0.15, 

0.25, and 0.35 were studied using the potential given by Eq.(30) with 
ifr « 0.06V and c ■ 10“ 4 mol/dm^. At this electrolyte concentration, the 
screening length <c is about 0.03 ym, and <d - 39.4. This results in the 
relatively short-ranged, repulsive potential shown in Figure 5. A nonaqueous 
suspension of sterically stabilized 1.2 u n diameter spheres was simulated at a 
volume fraction of 0.40 using the potential given by Eqs. (31) to (33). This 
potential has a deep primary minimum (not shown) for r/d < 1.005 and a shallow 
secondary minimum at r/d » 1.033 separated by a steeply repulsive wall that 
prevents irreversible coagulation into the primary minimum. The potential is 
shown in Fig. 5. 

Figures 6 through 9 contain the results for the aqueous suspensions. All 
calculations were performed using 32 particles and a time step of 2 x 10~ 4 s. 
Figure 6 compares the calculated g(r) at (> » 0.35 with the result of Van Megen 
and Snook.The agreement is very good. Figure 7 shows the progressive 
change in the equilibrium structure that occurs as the volume fraction 
decreases. The prominent nearest neighbor peak when $ - 0.35 arises from the 
strong repulsions that prevent particles from approaching each other too 
closely. This peak shrinks rapidly as $ decreases because the particles have 
more free volume in which to move. The excluded particle region around 
r/d - 1.55 is evidence that when $ * 0.35 the suspension has a high degree of 
solid-like order. This feature disappears as $ decreases. At $ » 0.25, g(r) 
is liquid-like, and at $ « 0.15, the suspension structure is beginning to 
resemble that of a repulsive soft-sphere gas. In Figs. 8 and 9, the m.s.d. 
and self-diffusion coefficient strongly show the effects of increasing $. Van 
Megen's and Snook's^ 0 results in Figs. 8 and 9 also illustrate the influence 
hydrodynamic interactions have in determining the short time diffusive motion 
of the particles. 
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Figure 6. Radial distribution function for dispersion with 
<d - 39.4 at 1 ■ 0.35 (from Ref. 30) . 
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Figure 



Plot ot the mean square displacement against time for dispersion 
with n - 0.1 mol Van Hagen's and Snooks' results with 

hydrodynamic interaction# (•«») and those neglecting 
hydrodynamic interaction* (-—) . The values alongside the curves 
denote the volume fraction. (After Ref. 10.) 



9. Plot of the time-dependent self-diffusion constant (Eg. (23)), 
expressed in units of D n , against time , other details as for 
Figure 8 (after Ref. 30). 
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Tha (laul4tlon of tha nonaquaoua auapanaion was dona uaing 108 particlaa 
and, dua to tha hardnaaa of tha rapulaiva potantial, a tiaa a tap of 5 * 10“ 6 a. 
In Fig. 10, g(r) for tha nonaquaoua auapanalon at $ - 0.4 la ahown. Tha 
aoat aignifleant faatura la tha vary atror.g naaraat naighbor paak. Thla paak 
la consldarably aharpar than thoaa found for tha aquaous auapanalona. Thla la 
dua to tha hardnaaa and ahortar ranga of tha nonaquaoua rapulaiva potantial In 
coablnatlon with tha aacondary attractiva wall. Tha paak thua rapraaanta 
waakly flocculatad particlaa with aaparatlona narrowly dlatrlbutad about tha 
aacondary alnlaun. Each partlcla haa, on avaraga, alx vary naar naighbor* 
with r/d lying batwaan 1.03 and 1.06. Tha aalf-diffuaion coafflclant for thia 
ayataa haa a vary rapid dacay (not ahown), dropping to 47 parcant of D 0 aftar 
only I0' 4 a and to 12 parcant of Dq aftar 1Q" 3 b. Tha asymptotic valua of about 
0.09 D 0 la raachad aftar 3 * 10* 3 a. Thla da cay la aaich faatar than that aaan 
In Fig. 9 for tha aquaoua auapanalona and la probably a conaaquanca of caging 
dua to tha waak flocculation or parhapa to tha alowar collactiva notion of 
waakly bound groupa of particlaa. 



r/d 


*-5343 


Figura 10. Equlllbrluai radial dlatributlon function for nonaquaoua 
dlaparalon with a » 0.4 
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4.2 NONEQUILIBRIUM SIMULATIONS 


Simula tions war# perforaed with 108 particlaa Cor voluat fractions of 0.2 
or 0.4, axcapt for ona run with + - 0.35, at various ahaar ratas ranging froa 
5 to 200s'" 1 . Tha particia diaaatar was always 1.2 un> Thraa diffsrant poten- 
tiais wara uaad in ordar to sssasa tha affact of intarparticia forcas on shaar 
thinning. Two of thasa wara also uaad for tha equilibrium siaulations just 
dascribad. Tha third potantial is givan by Eg. (30) with * - 0.06V and with 
c - IQ" 6 aol/da J . At this aisctroiyta concantration tha Oabya scraaning 
iangth is incraasad by a factor of tan ovar tha valua usad for tha equilibrium 
siaulations. With this incraasa, <d now equals 3.94, and tha ranga of tha 
rapuisiva, alactroatstic stabilisation potantial is substantially axtandad. 
This is apparant in Fig. 5, whars all thraa potantials ara plottad. 

Generally, tha first fiva to tan thousand tiaa stapa of a run wars usad 
to bring tha syataa to staady stats. Tha run would than ba contlnuad for 
another twenty to thirty thousand tiaa stapa, during which g(r), tha average 
of tha particia stress tensor (Eq. (26)) and its square wara coaputad. Tha 
coaputations ware perforaed on PSI'S in-house DEC MicroVax II coaputar system. 
In ordar to reduce tha computational tiaa needed for thasa preliainary stud¬ 
ies, randoa numbers with a unifora, rather than Gaussian, distribution wars 
generated over tha interval (-1/2, *1/2). A test run showed this to have only 
a small affact on tha calculatad results. For 108 particles, 1700 tiaa steps 
could ba perforaed par CPU hour. A typical 20,000 step run would than require 
about 12 hrs of coaputar tiaa. For tha nonaqueous systaaa, the same tiaa step 
(5 x 10" 6 s) uaad in tha equilibrium siaulations proved satisfactory. A time 
step of 2 x I0" 5 e was usad for siaulating tha aqueous suspensions. This is 
ona-tenth of tha valua usad for tha equilibrium siaulations. The saaller time 
step was needed to avoid artificially bringing pairs of particles into highly 
repulsive configurations that would result in unphysically large displacements 
on tha subsequent tiaa step. This can occur because with tha relatively soft 
rapuisiva potantials usad for tha aqueous systems studied hare, the effect of 
shaar is to drive soma of tha particles closer together where they experience 
stronger repulsions than at equilibrium. 
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This can be seen in Fig. 11 where the equilibrium g(r) and the 
angular-averaged nonequilibrium g(r) are plotted for an aqueous suspension 
with $ - 0.35 and <d - 39.4. Besides driving particles closer together, shear 
reduces the number of particles in the nearest neighbor shell, spreading them 
over a wider range of pair separations. This is evident from the appearance 
of particles in the range 1.45 < r/d < 1.65 where they were strongly excluded 
at equilibrium. Thus the average particle distribution is more uniform than 
at equilibrium, although it should be noted that in a sheared suspension g(r) 
is strongly angular dependent and not spherically symmetric as at equilibrium. 



Figure 11. Radial distribution functions for dispersion with < - 39.4 at 
t • 0.35, at equilibrium and at a shear rate of 10s~^ . 


The behavior of the sheared nonaqueous suspensions is somewhat different. 
Figures 12 through 14 show the angular-averaged nonequilibrium g(r) for 
) - 0.4 at three different shear rates: 20, 100, and 200s -1 . The very hard 
repulsive part of the potential for these systems prevents any significant 
decrease in the minimum pair separation at the shear rates studied. The high 
steric barrier also makes it virtually impossible for shear-induced coagulation 
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Figure 1 


Figure 1 



2. Nonequilibrium radial distribution function for aonaqueous 
dispersion with <*> ■ 0.4 and shear rate of 20s~^ . 



3. 


Nonequilibrium radial distribution function for nonaqueous 
dispersion with $ 0.4 and shear rate of 100s~^ . 
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Figure 14. Nonequilibrimn radial distribution function for nonaqueous 
dispersion with <j - 0.4 and shear rate of 200s~^ 

to occur at these shear rates. Thus the only significant effect of shear 
observed here is to reduce substantially the number of nearest neighbors 
around any particle. This is demonstrated by the marked reduction in height 
of the major peak as the shear rate increases. Compared to the equilibrium 
g(r) in Fig. 10, the sheared particle distribution is more uniform at larger r 
as well. Minor structural features in the equilibrium curve are steadily 
reduced with increasing shear rate. 

Figures 15 through 22 illustrate the calculated volume fraction and shear 
rate dependence of the relative viscosity of the various suspensions. Some 
comparisons with experimental data are also given. The relative viscosity n r 
of the suspension is defined as the ratio of the measured or calculated sus¬ 
pension viscosity to the viscosity of the pure suspending fluid, 

n r “ n/n 0 • (34) 
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Figure 15. Volume fraction dependence of relative viscosity for hard sphere 
suspensions (from Ref. 8) 
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Figure 16. The non-Newtonian behavior of nonaqueous latices at high volume 
fractions . Open symbols: data from Ref. 37. Solid symbols: 
nebd results for nonaqueous dispersion. 
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Figure 17. Relative viscosity versus Peclet number for sterically stabilized 
monodisperse polyvinyl chloride spheres in several solvents at 
it ° 0.20 (Ref. 38) . 
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Figure 20. Effect of added electrolyte on the relative viscosity of 
polystyrene latex suspension with a ■ 0.11 u°> at & - 0.40 
(Ref. 40): o deionized; □ 1.876 x 10~ 4 M HC1; 

'■1.876 x 10~-*M HC1; A 1 .876 ItT^M HC1; • 9.378 x 10~ 2 M 
HC1. (From Ref. 2). 

Figure 15 is adapted from the paper of Brady and Bossis. 8 It shows the volume 
fraction dependence of nr for several experimental systems, theoretical 
results of Beenakker, 41 and their Stokesian dynamics (SD) simulations. Our 
low shear (10s” 1 ) results for aqueous suspensions with <d - 39.4 at $ « 0.2, 
0.35, and 0.4 have been added to the figure. A higher shear rate (100s -1 ) 
value for $ - 0.4 is also included. The figure was originally intended to 
compare values for hard sphere suspensions, thus our values are not strictly 
comparable with the others. However, they do demonstrate that NEBD simula¬ 
tions of moderately repulsive spheres do produce a dependence on <j> that corre¬ 
sponds well with observed behavior on related systems. Our higher shear rate 
value for $ - 0.4 shows a considerable shear thinning effect; more discussion 
of this point will follow shortly. The lower line added to the figure repre¬ 
sents the relative viscosity due only to the solvent and independent particle 
hydrodynamic contributions. It neglects interparticle forces and hydrodynamic 
interactions. It is clear that interparticle forces contribute to our NEBD 
results to about the same degree that hydrodynamic interactions contribute to 
the force-free SD results of Brady and Bossis. 8 Brady and Bossis 8 also simu¬ 
lated a system with a much shorter-ranged, more highly repulsive potential 
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than used in any of the present aqueous simulations. This system showed a 
very weak shear-thinning effect, and exhibited a slight shear thickening at 
much higher shear rates due to hydrodynamic interactions. 

Figure 16 presents the volume fraction dependence of rir Cor nonaqueous 
latices in the high and low shear rate limit. 33 Also shown are our NEBD 
results for $ - 0.2 and 0.4 at shear rates of 20 and 200s** 1 . The discrepancy 
at ^ * 0.4 is almost entirely due to the lack of hydrodynamic interactions 
between particles in the simulations, but the shear thinning effect, due to 
interparticle forces, is reproduced reasonably well. Hydrodynamic interac¬ 
tions are relatively unimportant for the more dilute dispersion, which also 
shows very little shear thinning at these shear rates. Figure 17 is a compen¬ 
dium of data 38 on nonaqueous dispersions at $ » 0.2 taken from a review by 
Russel. 3 when plotted against a dimensionless shear rate, the Peclet number 
the relative viscosity data for different particle sizes lie on a 
single curve. Shear thinning is evident at very low dimensionless shear 
rates. This apparent conflict with the experimental results shown in Fig. 16 
may be due to some relatively weak, long ranged interaction or structure for¬ 
mation that is important only at very low shear rates. The available NEBD 
results plotted in Fig. 17 are in reasonable agreement with the data, but it 
was not possible to explore the shear thinning region. Meaningful results 
could not be obtained at lower shear rates due to large statistical fluctua¬ 
tions for the 108 particle system studied. Even for the two points shown, the 
standard deviation for the average particle stress contribution to n r was com¬ 
parable to the average value. 

Figure 18 shows the shear rate dependence of the two nonaqueous disper¬ 
sions calculated from the NEBD simulations. The error bars indicate the cal¬ 
culated standard deviations. It is gratifying that the NEBD simulations of 
the model nonaqueous system correspond well with data for real systems. It is 
also worth noting how much lower the viscosity of veil-stabilized dispersions 
is compared to that of carbon/hydrocarbon fuel slurries, such as reported by 
Lin and Brodkey. 38 An example of their results is shown in Fig. 19. These 
systems apparently have either a long-range network structure or a strong 
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shorter-ranged structure that leads to very high viscosities at low shear 
rates, to pronounced shear thinning behavior, and to high relative viscosities 
(«100) even at shear rates greater than 100s -1 . It seems reasonable to inter 
that improved polymeric stabilization techniques could greatly reduce the vis¬ 
cosities of these fuel slurries. 

Although it was not possible to simulate either long or short range 
structural effects in nonaqueous suspensions during the Phase I program, a 
crude analog for this, termed the secondary electroviscous effect, was simu¬ 
lated in the aqueous system studies. This effect refers to the increase in 
viscosity of aqueous dispersions that occurs when the electrolyte concentra¬ 
tion is progressively reduced. The phenomenon is illustrated with experi¬ 
mental results 4 ® in Fig. 20. At fixed particle volume fraction, the increas¬ 
ing Debye screening length results in much longer ranged forces, and the 
system becomes much more strongly interacting. Although this is clearly not 
the same type of strong interaction or long range structural effect assumed to 
be acting in the highly viscous nonaqueous slurry, it is nevertheless a valu¬ 
able confirmation of the applicability of the NEBD method to a strongly inter¬ 
acting system with long ranged forces. 

The results of the calculations are shown in Fig. 21. Aqueous disper¬ 
sions with <j> - 0.2 and <d = 39.4 and 3.94 were simulated. The change in < is 
due to a hundred-fold reduction in the salt concentration. The minimal shear 
thinning exhibited by the system with the much shorter ranged potential cor - 
trasts sharply with the behavior of the dispersion with long ranged forces. 

The shear rate dependence of n r for this latter system is so strong that below 
shear rates of 10s -1 , its viscosity exceeds that of a suspension at twice the 
volume fraction (<j> = 0.4) with the shorter ranged repulsive forces (<d » 

39.4). These latter results are plotted in Fig. 22. This system exhibits 
noticeable shear thinning at $ = 0.4, compared to the. minimal effect observed 
at l * 0.2, and the effect here is due to crowding particles more closely 
together. Reducing the interparticle spacing at fixed range of potential is 
roughly equivalent to increasing the range of the potential at fixed volume 
fraction. Both changes result in a more strongly interacting system, and that 
gives rise to an increase in the viscosity. 
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Based on the small number of cases studied, the following tentative con¬ 
clusions may be offered regarding the effect of interparticle forces on shear 
thinning. With increasing shear, the average number of near neighbors to any 
particle decreases and the average number of particles with larger separations 
grows. Both trends result in a decreased contribution to the average particle 
stress term and, thus, to a lower viscosity. Dispersions with harder, short 
ranged forces (even with long ranged, but weak attractions) show less shear 
thinning than ones with longer ranged strong interactions simply because the 
former systems have less "viscosity" to lose at any volume fraction. The 
system whose dominant interactions are short ranged generates significant con¬ 
tributions to the average particle stress only from very close pairs. In the 
system with a longer ranged strong interaction, the stress receives contribu¬ 
tions from a broader range of particle separations and, consequently, is 
larger. At sufficiently high shear rates, the particle distribution is so 
spread out that the force contribution to the stress is dominated by convec¬ 
tion. A complete discussion of shear thinning would include the effects of 
hydrodyamic interactions, 7 ’ 8 which become more important at higher shear 
rates. 44 Since these were not included in the Phase I feasibility study, no 
speculation will be offered regarding them. 

4.3 CLASSICAL STABILITY ESTIMATES 

Classical colloid stability estimates are based on the application of 
Fuchs' well known theory 45 for the rate of diffusion of particles over a 
potential barrier. The ratio of the rates without and with the barrier is 
termed the stability ratio, W. This is defined as 

W - d / exp(- i ^-)/r 2 dr (35) 

The nonaqueous system is of primary interest here, and ,the appropriate poten¬ 
tial is given by Eqs. (31) through (33). For parameter values of interest, 
this potential is very sharply peaked at a value of r, say r m , that is just 
slightly greater than d. Because of this, the exponential function in 
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Eq. (35) acts like a delta function, and a quadratic expansion of u(r) is suf¬ 
ficient to accurately approximate the integral via a Gaussian quadrature. The 
result is 

W - exp(u m /kT)/(d ✓(v" m /2irkT)) (36) 

where u m = u(r m ) and 

v" m - -(d 2 u/dr 2 )jr «r m 

Consider how stability is affected by the strength of the repulsive and 
attractive parts of the potential. For this exercise, the adsorbed polymer 
concentration C p will be used as the key parameter controlling the strength of 
the repulsive steric potential. The Hamaker constant A (Eq. (33)) is the nat¬ 
ural choice for the attractive potential. This quantity depends on the nature 
of two interacting particles as well as on the intervening fluid medium, and 
is not a free parameter. Tabor 46 indicates that metallic particles may have A 
values 5 to 10 times larger than those found for dielectric materials, so 
variations in A may be regarded as corresponding to changes in the material 
composition of the dispersed particles. 

Smoluchowski 1 s theory of rapid coagulation may be used to estimate the 
coagulation time t c , 

t c - 3n 0 /(4kTp) 

With the value p » 10^ 2 /cm® for 1 pm diameter particles at a volume fraction 
of 0.4, the coagulation time is of the order of 0.2s. This coagulation time 
must be increased by about a factor of 10® if stability for a year (3 x 10 7 s) 
is desired. Consequently, when W exceeds 10®, reasonable stability is 
insured; for W smaller than 10®, coagulation proceeds at a faster rate. 

Because p varies as d“3 a t fixed $, smaller particles require larger values of 
W to insure comparable periods of stability. It is well known that suspen¬ 
sions of small particles are generally more difficult to stabilize than those 
of larger particles. 22 
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Values of the stability ratio are shown in Table 1 for different combina¬ 


tions of Cp and A. Not surprisingly, the results show that decreases in Cp or 
increases in A result in decreased colloidal stability. Less adsorbed polymer 
means a less effective steric barrier. An increased attractive potential also 
reduces the repulsive barrier height. In a given system, A will be fixed, and 
stability can be enhanced by increasing the amount of polymer adsorbed. It 
should be noted that A depends, in principle, on the amount of adsorbed poly¬ 
mer, but this effect should be negligible for particle sizes of interest due 
} to the relatively small amounts of polymer involved. 

TABLE 1. Stability Ratio W, Potential Height u ra , and curvature v m " at r m 
for Different Concentrations of Adsorbed Polymer and Hamaker 
Constants. 


, Cp 3 

( 1 0~ 2 g/cm 3 ) 

A 

( 10 - 20 J) 

Un/kT 

d 2 v m 

'/2kT 

w 

Stable/ 

Unstable 

1 .3 

5.0 

154.2 

3.38 

X 

10 6 

9.0 

X 

10 63 

S* 

1 .0 

5.0 

61 .3 

1 .39 

X 

10 6 

6.4 

X 

1 0 23 

s 

0.95 

5.0 

49.4 

1.38 

X 

10 6 

4.3 

X 

10 18 

s 

0.90 

5.0 

38.5 

1.14 

X 

10 6 

8.7 

X 

1013 

s 

0.85 

5.0 

28.6 

0.83 

X 

10 6 

5.1 

X 

10 9 

s 

0.80 

5.0 

19.8 

0.71 

X 

10 6 

8.0 

X 

10 5 

u 

0.75 

5.0 

11.9 

0.61 

X 

10 6 

3.3 

X 

10 2 

u 

0.75 

1 .0 

64.5 

1 .71 

X 

10 6 

1 .4 

X 

1 0 2S 

s 

0.65 

1 .0 

43.1 

0.90 

X 

10 6 

3.2 

X 

10 14 

s 

0.55 

1.0 

25.8 

0.53 

X 

10 6 

4.0 

X 

10 8 

s 

0.45 

1 .0 

12.6 

0.33 

X 

10 6 

8.9 

X 

10 2 

u 

0.40 

1 .0 

7.5 

0.20 

X 

10 6 

7.0 



u 

0.40 

0.5 

14.2 

0.34 

X 

10 6 

4.3 

X 

10 3 

u 

0.50 

0.5 

27.4 

0.61 

X 

10 6 

1.9 

X 

10 9 

s 

0.55 

0.5 

35.6 

0.86 

X 

10 6 

5.5 

X 

1012 

s 


* System used for computer simulations. 
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5 . 


SUMMARY AND CONCLUSIONS 



A number of equilibrium and nonequilibrium Brownian dynamics (NEBD) simu¬ 
lations of concentrated colloidal suspensions have been performed. Both 
aqueous and nonaqueous suspensions have been studied. Stability estimates for 
sterically stabilized nonaqueous suspensions have been made. These results 
can be understood in terms of the changes occurring in the repulsive steric 
potential and attractive van der Waals potential as the adsorbed polymer con¬ 
centration and Hamaker constant are changed. Quantitative evidence for the 
disturbance of the equilibrium structure of the dispersions by shearing has 
been obtained. Shear viscosities of suspensions have been calculated as a 
function of shear rate and particle volume fraction. With allowance for the 
neglect of hydrodynamic interactions, good agreement was obtained with experi¬ 
mental viscosities for comparable systems. The simulation results are sensi¬ 
tive to the type, strength and range of the potential interactions used to 
describe the dispersion. The overall conclusion of this Phase I study is that 
the NEBD technique is capable of describing the essential features of sheared 
suspension behavior. These Phase I results establish the feasibility and des¬ 
irability of using NEBD simulation methods as a predictive engineering tool 
for the design of slurries. 


40 




6. REFERENCES 


1. Jeffrey, D.J. and Acrivos, A., AIChE Journal _22, 417 (1976). 

2. Russel, W.B., in Theory of Dispersed Multiphase Flow (Academic, London 
1983) , p. 1 . 

3. Goodwin, J.W., in Colloid Science, Volume 2 (Specialist Periodical 
Reports) ed. D.H. Everett (The Chemical Society, London, 1975), p246. 

4. Goodwin, J.W. , in Colloidal Dispersions , ed. J.W. Goodwin (Royal Society 
of Chemistry, London 1982), p. 165. 

5. Hirtzel, C.S. and Rajagopalan, R., Colloidal Phenomena (Noyes, Park 
Ridge, New Jersey) 1985). 

6. Ansel, G.C., Dickinson, E., and Ludvigsen, M., J. Chem. Soc., Faraday 

Trans. 2, j3J_, 1269 (1985). 

7. Bossis, G. and Brady, J.F., J. Chem. Phys. _80, 5141 (1984). 

8. Brady, J. and Bossis, G., J. Fluid Mech., 155 , 105 (1985). 

9. Hoover, W.G., Ann. Rev. Phys. Chem. _34, 103 (1983). 

10. Evans, D.J., Molec. Phys. 42 , 1355 (1981); _37, 1745 (1979). 

11. Evans, D.J. and Morriss, G.P., Comput. Phys. Repts. J., 299 (1984). 

12. Heyes, D.M., Montrose, C.J., and Litovitz, T.A., J. Chem. Soc. Faraday 
Trans. 2, 79, 61 1 ( 1983) . 

13. Hoover, W.G., Molecular Dynamics , Lecture Notes in Physics 258 (Springer- 
Verlag, Berlin, 1986). 

14. Ermak, D.L., J. Chem. Phys. _62, 4189, 4197 (1975). 

15. Doll, J.D. and Dion, D.R., Chem. Phys. Lett. _37, 386 (1976). 

16. Turq, P., Lantelme, F., and Friedman, H.L., J. Chem. Phys. _66, 3039 
(1977). 

17. Ermak, D.L. and McCammon, J.A., J. Chem. Phys. _69, 1352 (1978). 

18. Fixman, M., J. Chem. Phys. _69, 1527, 1538 ( 1978); Macromolecules _1_4, 1710 
(1981) . 

19. Israelachvili, J.N., Adv. Coll. Interface Sci. 16, 31 (1982). 


41 








20. Ninham, B.W., Adv. Coll. Interface Sci. _1_6, 3 (1982). 

21. Napper, D.H., Colloidal Dispersions , ed. J.W. Goodwin (Royal Society of 
Chemistry, London, 1982), p. 99. 

22. Verwey, E.J.W. and Overbeek, J. Th. G., The Theory of the Stability of 
Lyophobic Colloids (Elsevier, Amsterdam, 1948). 

23. Sonntag, H. and Strenge, K., Coagulation and Stability of Disperse 
Systems (Israel Program for Scientific Translations, Hlasted Press, 

New York, 1972). 

24. Gregory, J., J. Coll. Interface, Sci. _£3, 138 (1981). 

25. Wilemski, G. , J. Coll. Interface, Sci. 88 _, 111 (1982). 

26. Wilemski, G., J. Stat. Phys. _1_4, 153 ( 1976). 

27. Fincham, D., Comput. Phys. Commun. _21_, 247 (1980). 

28. Lees, A. and Edwards, S., J. Phys. C_5, 1921 (1972). 

29. Van Megen, W. and Snook, I., Faraday Discuss. Chem. Soc. _76, 151 (1983). 

30. Van Megen, W. and Snook, I., J. Chem. Soc., Faraday Trans. 2, _80, 383 

(1984). 

31. Cook, R., Private communication (1986). 

32. Ball, R.C. and Richmond, P. , Phys. Chem. Liq. 99 (1980). 

33. Batchelor, G.K., J. Fluid Mech. 83 ., 97 (1977). 

34. Ottewill, R.H. and Walker, T., Kolloid Z.Z. Polym. 227 , 108 (1968). 

35. Cairns, R.J.R., Van Megan, W., and Ottewill, R.H., J. Colloid Interface 
Sci. 2i. 511 (1981) . 

36. Hamaker, H.C. Physica 1058 (1937). 

37. Papier, Y.S. and Krieger, I.M., J. Colloid Interface Sci. _34, 126 (1970). 

38. Willey, S.J. and Macosko, C.W., J. Rheol. _22, 525 (1978). 

39. Lin, S.-F. and Brodkey, R.S., J. Rheol. 29 ., *47 (1985). 

40. Krieger, I.M. and Equiluz, M., Trans. Soc. Rheol. _20, 29 (1976). 

41. Beenakker, C.W.J., Physica A 128 , 44 (1984). 


42 





42. (a) Gaylor, K.J., Snook, I.K., Van Megen W.J., and Watts, R.O., Chern. 
Phys. 43, 233 (1979); (b) J.C.S. Faraday II 76. 1067 (1980); (c) J. Phys. 
A J_3, 2513 ( 1980) . 

43. Snook, I.K., Van Megen, W., Gaylor, K.J., and Watts, R.O., Adv. Coll. 
Interface Sci. 17 , 33 (1982). 

44. De Kruif, C.G., Van Iersel, Vrij, A., and Russel, W.B., J. Chem. Phys. 

83 , 4717 (1985). 

45. Fuchs, N.A., Z. Phys. Chem. A171 , 199 (1934). 

46. Tabor, D., in Colloidal Dispersions , ed. J.W. Goodwin (Royal Society of 
Chemistry, London 1982), p. 23. 


43/44 


APPENDIX A 


EQUILIBRIUM SIMULATIONS OF DILUTE DISPERSIONS 

Gaylor, Snook, Van Megen and Watts 2 ^ ■30,42,43 have published the results 
of a number of equilibrium Brownian dynamics simulations of colloidal suspen¬ 
sions. In order to verify the satisfactory operation of the present equili¬ 
brium code, several simulations were performed for conditions reported in the 
literature. The potentials used are shown in Figure A-1. Comparisons were 
made for the radial distribution function, mean square displacement and self¬ 
diffusion coefficient. Dilute, but strongly interacting, suspensions were 
studied. Hydrodynamic interactions were neglected. This is justifiable for 
the dilute systems. 

Comparisons of the present computations for dilute suspensions with those 
of Gaylor et al. 423 ' 4 ^ anc j Van M e gen and Snook 2 ^ are shown in Figs. A-2 
through A-6. Generally the calculations were done using 32 particles; however 
those at $ = 9.9 x 10 -4 (Figs. A-5 and A-11) required 256 particles. In 
general the comparisons are quite favorable. The results of Gaylor et al. 42a 
in Fig. A-2 show how changing the ionic strength of the electrolyte affects 
the structure of the suspension. At the highest concentration, the repulsive 
potential is relatively short-ranged compared to the average interparticle 
spacing, and g(r) is characteristic of a very dilute gas with a soft, purely 
repulsive potential. As the electrolyte concentration decreases, the range of 
the potential increases. Because the particle volume fraction is constant, 
the suspension behaves as though it were becoming denser, and g(r) develops 
its characteristic nearest-neighbor peak due to the excluded volume around 
each particle. In other words, particles are being pushed away from each other 
as the repulsive force between them increases. The effect of volume fraction 
on the time dependence of the m.s.d. and the self-diffusion coefficient as 
calculated by Gaylor et al. is shown in Figs. 8 and 9. As expected, motion 
slows with increasing $. These two figures show clearly how the average 
motion differs at short and long times. Because the average force on a parti¬ 
cle is zero, short time motion will be force-free diffusion as described by 
Eq. (19) with a self-diffusion coefficient equal to D^. Over long times, 
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Figure A-1. Repulsive potentials used for simulations of dilute dispersions . 

Potentials are based on Eq. (29) and are due to electrical double 
layer interactions in an electrolyte with a concentration of 
10 -6 mol/dm 3 . Surface potential i^ 0 is 0.22V when d = 0.046 ym 
and 0.23V when d = 0.09 yin. 



Figure A-2. Radial distribution function as a function of electrolyte concen ¬ 
tration when <{) = 1.5 x 10 -4 ; n ~ 10 -3 mol m” 3 , — • — n = 

10”^ mol ra“ J , - n » 10 - ^ mol ra -3 (after Ref. 42a). 
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because of the slowly changing particle configuration, forces on the particles 
are correlated, and the average motion is strongly affected. 

Our test calculations in Figs. A-2 through A-4 were made for a volume 
fraction of 1.5 x 10 -4 using two different fluid viscosities, 1 and 2 mPa*s. 
Doubling the usual value of ti 0 results in halving the single particle diffu¬ 
sion coefficient (cf. Eq. (11)). This significantly retards the average 
motion of particles in suspension, as can be seen in Figure A-3. In princi¬ 
ple, changing no should not affect g(r) because it is an equilibrium property 
solely determined by U. The results in Fig. A-2 approximately bear this out. 
The discrepancy may be due to using too large a time step At in the simula¬ 
tions. Van Megen and Snook 30 have noted that if At is too large, peaks in 
g(r) tend to be artificially flattened. Doubling the viscosity while holding 
At constant effectively halves the time step, and the peak height grows as a 
consequence. It should be noted that the same At reported by Gaylor et al. 42a 
was used in these calculations. The remaining discrepancy with their results 
may be due to running the simulations for an insufficient number of time 
steps. 

In Figs. A-5 and A-6, the results of Van Megen and Snook 29 show the 
effect of increasing the volume fraction to the value 9.9 x 10 -4 . Their 
results in Fig. A-5 using 256 particles show that the suspension is becoming 
increasingly ordered with the particles distributed about the sites of a bcc 
lattice. The strong interactions that promote ordering also greatly inhibit 
particle motion. This can be seen in Fig. A-6 from the slowly increasing 
m.s.d. at the highest volume fraction. Our results for g(r) agree well at 
r/d<14, but are poorer for larger values of r. This is probably due to our 
neglect of long range contributions to the total force from particles beyond 
the cut-off distance of 1 / 2 . The calculated m.s.d. agrees well probably 
because its early behavior should be controlled by interactions with particles 
in the first and second nearest neighbor shells. 





